index_upload

1. Data preparation

library(dplyr)
library(traitstrap)
library(ggplot2)

morphometric_data <- utils::read.csv(here::here("data", "morphometric_data.csv"), sep = ";", header = T, dec = ".")

morpho_data <- morphometric_data %>%
  select(-c(variable_abbreviation, variable_unit)) %>%
  t() %>%
  as.data.frame() %>%
  janitor::row_to_names(row_number = 1) %>%
  `rownames<-`(NULL)%>%
  # delete for now (n=1) and Eurypharynx pelecanoides present 7 missing traits 
  filter(!species%in% c("Diaphus_sp","Eurypharynx_pelecanoides"))

# replace empty value by NA 
morpho_data[morpho_data ==""] <- NA

# Numeric variables
morpho_data[, 4:23] <- sapply(morpho_data[, 4:23], as.numeric)

1.1.Data imputation

• mice algorithm: n imputation = 5, n iterations = 50

#select numeric variables for imputation 
original_data <- morpho_data %>%
   filter(species!="Eurypharynx_pelecanoides") %>% 
  select(1:23)

imputation <-
  mice::mice(
    original_data,
    m = 5,
    maxit = 50,
    printFlag = F
  )

imputed_data <- mice::complete(imputation)

1.2. Species * traits

• calculate functional traits

# calculate functional numeric traits
numeric_traits <- imputed_data %>%
  na.omit() %>%
  select(-individual_code) %>%
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume
  ) %>%
  group_by(species) %>%
  summarise(across(everything(), mean, na.rm = TRUE)) %>%
  arrange(species)

# categorical traits for species without NA
cat_morpho <- morpho_data %>%
  select(
    species,
    ventral_photophores,
    gland_head,
    chin_barbel,
    small_teeth,
    large_teeth,
    fang_teeth,
    retractable_teeth,
    internal_teeth,
    gill_raker_types,
    oral_gape_axis
  ) %>%
    na.omit() %>%
  distinct() %>%
  arrange(species)

# combined the two data frames
fish_traits <- numeric_traits %>%
  inner_join(cat_morpho, by = "species") %>%
  arrange(species) %>% 
  mutate(species = stringr::str_replace(species,  "Cyclothone_sp","Cyclothone")) %>% 
  filter(species!="Eurypharynx_pelecanoides") %>% 
  tibble::column_to_rownames("species")%>%
  # assign trait type 
  # as.factor for qualitative traits
  mutate_if(is.character, as.factor)%>%
  # as.ordered for ordered variables
  mutate_at(c("gill_raker_types", "oral_gape_axis"), as.ordered)

Data summary

fish_traits_sum <- fish_traits %>% 
  mutate(across(where(is.numeric), round, 2))

## Display the table ----
htmltools::tagList(DT::datatable(fish_traits_sum))

traits correlation

M <-cor(numeric_traits[, c(-1)])

ggcorrplot::ggcorrplot(M, hc.order = TRUE, type = "lower",
                       lab = TRUE, tl.cex = 9, lab_size = 3)

# list of species 
sp_names <- c(rownames(fish_traits))

# taxonomic_families
taxonomic_families <- sp_names %>%
  as.data.frame() %>%
  `colnames<-`("species") %>% 
  mutate(
    family = case_when(
      species %in%
        c(
          "Benthosema_glaciale",
          "Ceratoscopelus_maderensis",
          "Diaphus_metopoclampus",
          "Lampanyctus_ater",
          "Lampanyctus_crocodilus",
          "Lampanyctus_macdonaldi",
          "Lobianchia_gemellarii",
          "Myctophum_punctatum",
          "Notoscopelus_bolini",
          "Notoscopelus_kroyeri",
          "Bolinichthys_supralateralis"
        ) ~ "Myctophidae",
      species %in% c(
        "Borostomias_antarcticus",
        "Chauliodus_sloani",
        "Malacosteus_niger",
        "Melanostomias_bartonbeani",
        "Stomias_boa"
      ) ~ "Stomiidae",
      species %in% c(
        "Holtbyrnia_anomala",
        "Holtbyrnia_macrops",
        "Maulisia_argipalla",
        "Maulisia_mauli",
        "Maulisia_microlepis",
        "Normichthys_operosus",
        "Searsia_koefoedi",
        "Sagamichthys_schnakenbecki"
      ) ~ "Platytroctidae",
      species %in% c("Sigmops_bathyphilus",
                     "Gonostoma_elongatum") ~ "Gonostomatidae",
      species %in% c(
        "Argyropelecus_hemigymnus",
        "Maurolicus_muelleri",
        "Argyropelecus_olfersii"
      ) ~ "Sternoptychidae",
      species == "Anoplogaster_cornuta" ~ "Anoplogastridae",
      species %in% c("Arctozenus_risso", "Paralepis_coregonoides") ~ "Paralepididae",
      species == "Bathylagus_euryops" ~ "Bathylagidae",
      species == "Cyclothone" ~ "Gonostomatidae",
      species == "Derichthys_serpentinus" ~ "Derichthyidae",
      species == "Eurypharynx_pelecanoides" ~ "Eurypharyngidae",
      species == "Evermannella_balbo" ~ "Evermannellidae",
      species == "Lestidiops_sphyrenoides" ~ "Lestidiidae",
      species == "Melanostigma_atlanticum" ~ "Zoarcidae",
      species %in% c("Photostylus_pycnopterus",
                     "Xenodermichthys_copei") ~ "Alepocephalidae",
      species == "Serrivomer_beanii" ~ "Serrivomeridae"
    )
  )%>% 
  mutate(
    order = case_when(
      family =="Myctophidae" ~ "Myctophiformes",
      family %in% c("Stomiidae","Gonostomatidae", "Sternoptychidae") ~  "Stomiiformes",
      family %in% c("Platytroctidae","Alepocephalidae") ~ "Alepocephaliformes",
      family == "Anoplogastridae" ~ "Trachichthyiformes",
      family %in% c("Paralepididae","Evermannellidae","Lestidiidae") ~ "Aulopiformes",
      family ==  "Bathylagidae" ~ "Argentiniformes",
      family %in% c("Derichthyidae","Serrivomeridae") ~ "Anguilliformes",
      family ==  "Eurypharyngidae" ~"Saccopharyngiformes",
      family == "Zoarcidae" ~ "Perciformes",
    )
  )

1.3. Species * assemblages matrix

Number of trawl hauls per depth

• Epipelagic = 8
• Upper mesopelagic = 26
• Lower mesopelagic = 16
• Bathypelagic = 16

# Metadata
metadata <-  utils::read.csv(here::here("data", "metadata.csv"), sep = ";", header = T, dec = ".")%>%
  # calculation of standardized biomass values (vertical  trawl opening * horizontal trawl opening * distance traveled)  
  mutate(volume_filtered = 24*58*distance_trawled_m)

ggplot(metadata, aes(x=depth))  +
  ylab ("Number of trawls")+
  xlab("Immersion depth (m)")+
  geom_histogram(binwidth=100, col="white", fill=alpha("black",0.55))+
  theme_light()+
  coord_flip()+ 
  scale_x_reverse()+
  labs(fill= "")+
  guides(fill="none")+
  scale_y_continuous(breaks = c(2,4,6,8,10,12))+
  theme(axis.text.x= element_text(size=12),
        axis.text.y= element_text(size=12),
        axis.title.y = element_text( size=12),
        axis.ticks = element_blank())

Biomass data

biomass_sum <- depth_fish_biomass %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "depth_layer") %>%
  tidyr::pivot_longer(!depth_layer, names_to = "species", values_to = "biomass") %>%
  group_by(depth_layer) %>%
  filter(biomass > 0) %>%
  summarise(biomass_depth = round(sum(biomass), 3),
            n = n())

htmltools::tagList(DT::datatable(biomass_sum))

1.4. Traits types

The first column contains traits name. The second column contains traits type following this code:

• N: nominal trait (factor variable)
• O: ordinal traits (ordered variable)
• Q: quantitative traits (numeric values)

fish_traits_cat <- utils::read.csv(here::here("data", "fish_traits_cat.csv"), sep = ";", header = T, dec = ".") 
htmltools::tagList(DT::datatable(fish_traits_cat))

2.Functional spaces

2.1 Compute data summaries

## Summary of the assemblages * species data.frame ----
asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)
asb_sp_fish_occ  <- asb_sp_fish_summ$"asb_sp_occ"

htmltools::tagList(DT::datatable(asb_sp_fish_occ))

2.2 Computing distances between species based on functional traits

• We have non-continuous traits so we use the Gower distance (metric = “gower”) as this method allows traits weighting.
• scale_euclid = TRUE

sp_dist_fish <- mFD::funct.dist(
  sp_tr         = fish_traits,
  tr_cat        = fish_traits_cat,
  metric        = "gower",
  scale_euclid  = "scale_center",
  ordinal_var   = "classic",
  weight_type   = "equal",
  stop_if_NA    = TRUE)

2.3 Building functional spaces and chosing the best one

2.3.1 Computing several multimensional functional spaces and assessing their quality

• mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).

fspaces_quality_fish <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 10,
  deviation_weighting = "absolute",
  fdist_scaling       = FALSE,
  fdendro             = "average")

## Quality metrics of functional spaces ----
round(fspaces_quality_fish$"quality_fspaces", 3)

               mad
pcoa_1d      0.143
pcoa_2d      0.079
pcoa_3d      0.050
pcoa_4d      0.030
pcoa_5d      0.023
pcoa_6d      0.017
pcoa_7d      0.014
pcoa_8d      0.015
pcoa_9d      0.017
pcoa_10d     0.019
tree_average 0.040

Variance explained by each axis

# Extract eigenvalues information
eigenvalues_info <- fspaces_quality_fish$"details_fspaces"$"pc_eigenvalues"

# Create a dataframe to store the results
variance_df <- data.frame(
  PC = c("PC1", "PC2", "PC3", "PC4"),
  VarianceExplained = c(
    eigenvalues_info[1, "Cum_corr_eig"] * 100,
    (eigenvalues_info[2, "Cum_corr_eig"] - eigenvalues_info[1, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[3, "Cum_corr_eig"] - eigenvalues_info[2, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[4, "Cum_corr_eig"] - eigenvalues_info[3, "Cum_corr_eig"]) * 100
  )
) %>% 
  mutate(across(where(is.numeric), round, 2))

htmltools::tagList(DT::datatable(variance_df))

2.3.2 Illustrating the quality of the functional spaces

mFD::quality.fspaces.plot(
  fspaces_quality            = fspaces_quality_fish,
  quality_metric             = "mad",
  fspaces_plot               = c("tree_average", "pcoa_2d", "pcoa_3d", 
                                 "pcoa_4d", "pcoa_5d", "pcoa_6d"),
  name_file                  = NULL,
  range_dist                 = NULL,
  range_dev                  = NULL,
  range_qdev                 = NULL,
  gradient_deviation         = c(neg = "darkblue", nul = "grey80", pos = "darkred"),
  gradient_deviation_quality = c(low = "yellow", high = "red"),
  x_lab                      = "Trait-based distance")

2.3.3 Testing the correlation between functional axes and traits

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

# As we have 26 traits we have to split the df to see correlation between functional axes and traits 
# first set ----
fish_traits_1 <- fish_traits%>%
  select(1:9)

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_1, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value" < 0.05), ]

                trait axis         test stat value p.value
1            eye_size  PC1 Linear Model   r2 0.164  0.0087
2            eye_size  PC2 Linear Model   r2 0.174  0.0067
3            eye_size  PC3 Linear Model   r2 0.101  0.0425
4            eye_size  PC4 Linear Model   r2 0.162  0.0091
5      orbital_length  PC1 Linear Model   r2 0.406  0.0000
6      orbital_length  PC2 Linear Model   r2 0.104  0.0394
10       gill_outflow  PC2 Linear Model   r2 0.513  0.0000
13  oral_gape_surface  PC1 Linear Model   r2 0.126  0.0225
14  oral_gape_surface  PC2 Linear Model   r2 0.500  0.0000
18    oral_gape_shape  PC2 Linear Model   r2 0.150  0.0123
24 oral_gape_position  PC4 Linear Model   r2 0.228  0.0016
26   lower_jaw_length  PC2 Linear Model   r2 0.685  0.0000
29        head_length  PC1 Linear Model   r2 0.330  0.0001
30        head_length  PC2 Linear Model   r2 0.382  0.0000
34         body_depth  PC2 Linear Model   r2 0.353  0.0000
35         body_depth  PC3 Linear Model   r2 0.196  0.0038

## Plot ----
fish_tr_faxes$"tr_faxes_plot"

# second set ----
fish_traits_2 <- fish_traits%>%
  select(10:18)

fish_tr_faxes_2 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_2, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value" < 0.05), ]

                    trait axis           test stat value p.value
2   pectoral_fin_position  PC2   Linear Model   r2 0.278  0.0004
5  pectoral_fin_insertion  PC1   Linear Model   r2 0.282  0.0004
6  pectoral_fin_insertion  PC2   Linear Model   r2 0.406  0.0000
9       transversal_shape  PC1   Linear Model   r2 0.116  0.0297
11      transversal_shape  PC3   Linear Model   r2 0.218  0.0021
14  caudal_throttle_width  PC2   Linear Model   r2 0.401  0.0000
19   dorsal_fin_insertion  PC3   Linear Model   r2 0.190  0.0044
23           eye_position  PC3   Linear Model   r2 0.191  0.0043
26       operculum_volume  PC2   Linear Model   r2 0.104  0.0393
32    ventral_photophores  PC4 Kruskal-Wallis eta2 0.590  0.0000
33             gland_head  PC1 Kruskal-Wallis eta2 0.245  0.0011

## Plot ----
fish_tr_faxes_2$"tr_faxes_plot"

# third set ----
fish_traits_3 <- fish_traits%>%
  select(19:26)

fish_tr_faxes_3 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_3, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value" < 0.05), ]

              trait axis           test stat value p.value
1       chin_barbel  PC1 Kruskal-Wallis eta2 0.142  0.0107
3       chin_barbel  PC3 Kruskal-Wallis eta2 0.116  0.0186
4       chin_barbel  PC4 Kruskal-Wallis eta2 0.176  0.0051
5       small_teeth  PC1 Kruskal-Wallis eta2 0.323  0.0002
9       large_teeth  PC1 Kruskal-Wallis eta2 0.454  0.0000
10      large_teeth  PC2 Kruskal-Wallis eta2 0.214  0.0022
13       fang_teeth  PC1 Kruskal-Wallis eta2 0.410  0.0000
21   internal_teeth  PC1 Kruskal-Wallis eta2 0.076  0.0467
23   internal_teeth  PC3 Kruskal-Wallis eta2 0.634  0.0000
25 gill_raker_types  PC1 Kruskal-Wallis eta2 0.331  0.0007
26 gill_raker_types  PC2 Kruskal-Wallis eta2 0.366  0.0004
30   oral_gape_axis  PC2 Kruskal-Wallis eta2 0.378  0.0003
31   oral_gape_axis  PC3 Kruskal-Wallis eta2 0.243  0.0036

## Plot ----
fish_tr_faxes_3$"tr_faxes_plot"

Summary of traits with a significant effect

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = F)

## Print traits with significant effect ----
traits_effect <- fish_tr_faxes[which(fish_tr_faxes$p.value< 0.05),] %>% 
  as.data.frame() %>% 
  arrange(axis, desc(value))

htmltools::tagList(DT::datatable(traits_effect))

2.4 Plotting the selected functional space and position of species

# Convert sp_faxes_coord_fish to a data frame
sp_coord_community <- as.data.frame(sp_faxes_coord_fish[, c("PC1", "PC2")]) %>%
  tibble::rownames_to_column(var = "species")

# Depth_layer ----
# Transform biomass data
sp_all_layers <- depth_fish_biomass %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "depth_layer") %>%
  tidyr::pivot_longer(!depth_layer, values_to = "biomass", names_to = "species") %>%
  filter(biomass > 0)

# Create a list to store data frames for each depth layer
all_data_layers <- lapply(unique(sp_all_layers$depth_layer), function(layer) {
  sp_layer <- sp_all_layers %>%
    filter(depth_layer == layer)
  
  # Mark presence or absence
  sp_coord_layer <- sp_coord_community %>%
    mutate(layer_presence = ifelse(species %in% sp_layer$species, "yes", "no")) %>%
    mutate(depth_layer = layer)
  
  return(sp_coord_layer)
})

# Combine all layers into one data frame
sp_coord_community_all_layers <- bind_rows(all_data_layers)

# Function to calculate the convex hull for a given data frame
calculate_hull <- function(data) {
  if (nrow(data) < 3) return(data)  # Handle cases with fewer than 3 points
  data %>%
    slice(chull(PC1, PC2))
}

# Initialize lists to store hulls
hull_all_combined <- NULL
hull_layer_combined <- NULL

for(layer in unique(sp_coord_community_all_layers$depth_layer)) {
  sp_layer <- sp_coord_community_all_layers %>%
    filter(depth_layer == layer)
  
  hull_all <- calculate_hull(sp_layer) %>%
    mutate(depth_layer = layer, type = "all")
  
  hull_layer <- calculate_hull(sp_layer %>% filter(layer_presence == "yes")) %>%
    mutate(depth_layer = layer, type = "present")
  
  hull_all_combined <- bind_rows(hull_all_combined, hull_all)
  hull_layer_combined <- bind_rows(hull_layer_combined, hull_layer)
}

# Define depth layer levels in desired order
depth_levels <- c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic")

# Ensure depth_layer is a factor with specified levels
sp_coord_community_all_layers$depth_layer <- factor(sp_coord_community_all_layers$depth_layer, levels = depth_levels)
hull_all_combined$depth_layer <- factor(hull_all_combined$depth_layer, levels = depth_levels)
hull_layer_combined$depth_layer <- factor(hull_layer_combined$depth_layer, levels = depth_levels)

# Define colors for each depth layer
depth_colors <- c("Epipelagic" = "#FEA520", 
                  "Upper mesopelagic" = "#D62246", 
                  "Lower mesopelagic" = "#6255B4", 
                  "Bathypelagic" = "#3C685A")

# Plotting
plot_depth <- ggplot() +
  # Plot present species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "yes"), 
             aes(PC1, PC2, col = depth_layer), size = 2, shape=19) +
  # Plot absent species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "no"), 
             aes(PC1, PC2), col = "grey", shape = 8, size = 2) +
  # Plot hulls
  geom_polygon(data = hull_all_combined, aes(x = PC1, y = PC2), 
               fill =  "gray70", alpha = 0.2) +
  geom_polygon(data = hull_layer_combined, aes(x = PC1, y = PC2, group = depth_layer, fill = depth_layer), 
               alpha = 0.4) +
  theme_light() +
  scale_color_manual(values = depth_colors) +
  scale_fill_manual(values = depth_colors) +
  guides(col = "none", shape = "none", fill = "none") +
  facet_wrap(~ depth_layer, nrow=1 ) +  
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )


beta_fd_indices_fish <- mFD::beta.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_occ       = asb_sp_fish_occ,
  check_input      = TRUE,
  beta_family      = c("Jaccard"),
  details_returned = T)

Serial computing of convex hulls shaping assemblages with conv1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_geom_coord

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%

vertices_community <-beta_fd_indices_fish$"details"$"asb_vertices"$Bathypelagic 

sp_coord_community <- as.data.frame(sp_faxes_coord_fish[, 1:2]) %>% 
  tibble::rownames_to_column(var = "species") %>% 
  mutate(vertices=case_when(species%in%vertices_community~"vertice",
                            !species%in%vertices_community~"not_vertice")) 

hull <- sp_coord_community%>% 
  slice(chull(PC1, PC2))

plot_com <- ggplot(data = sp_coord_community,aes(PC1, PC2)) +
  scale_color_manual(values = c("grey", "black"))+
  geom_point(size = 2.5, alpha=0.8, aes( shape=vertices, col= vertices)) +
  geom_polygon(data=hull, aes(x = PC1, y = PC2),
               fill = "gray70", alpha = 0.2) +
  theme_light() +
  geom_hline(yintercept = 0, linetype = 2, col = "gray45") +
  geom_vline(xintercept = 0, linetype = 2, col = "gray45") +
  guides(col = "none", shape = "none", fill = "none") +
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )

ggpubr::ggarrange(
  labels = c("A","B"),
  plot_com,   
  plot_depth, 
  nrow = 2
) 

ggsave("functional_space.png", path = "figures", dpi = 800, height = 8, width = 10)

3. SES Functional diversity indices

3.1.calculate and plot SES

• Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units

• null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurrence frequency)

• FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.

• FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.

• FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.

• FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.

• calcul des indices avec package FD

• correction distance matrix (obtenue avec Gower) = “none” or “lingoes” transformation avec “sqrt” ne permet pas de transformer la matrice en “euclidean”

• randomisation des traits des espèces par couche de profondeur

• interprétation : La richesse fonctionnelle (FRic) augmente avec la profondeur pour atteindre un maximum en bathypélagique, ce qui correspond à nos anciens résultats. Cependant, les valeurs de SES FRic sont également positives en épipélagique, ce qui indique une richesse fonctionnelle plus élevée qu’attendue au hasard. En revanche, les SES FDIS sont significativement plus faibles en épipélagique, ce qui indique que les espèces représentant le plus de biomasse sont similaires sur le plan fonctionnel.

• comment faire pour tester la symétrie et la normalité des valeurs de SES, suelement sur les quatres valeurs que j’obtiens par indices ?

# Model parameters ----
# Number of simulations
n_simulations <- 9

# correction method to use when the distance matrix cannot be represented in a Euclidean space
corr_method <- "lingoes"# "none" 

# Depth layers
depth_layers <- rownames(depth_fish_biomass)

# List of indices to be calculated
indices <- c("FRic", "FDis", "FDiv", "FEve") 

# Matrices for storing observed and simulated results
dbFD_result_obs <- matrix(NA, nrow = length(depth_layers), ncol = length(indices), 
                          dimnames = list(depth_layers, indices))

dbFD_result_sim <- array(NA, dim = c(length(depth_layers), n_simulations, length(indices)), 
                         dimnames = list(depth_layers, paste0("Sim.", 1:n_simulations), indices))

# Function to randomize traits while preserving factor levels----
randomize_traits <- function(traits_mat) {
  randomized_mat <- traits_mat
  # Randomisation of columns
  for (trait in colnames(traits_mat)) {
    randomized_mat[, trait] <- sample(traits_mat[, trait], replace = FALSE)
  }
  return(randomized_mat)
}

# Loop on each depth layer ----
for (layer in depth_layers) {
  
  # Select the species present in the selected depth layer
  species_in_layer <- colnames(depth_fish_biomass)[depth_fish_biomass[layer, ] > 0]
  
  # Filter lines (traits) and remove constant columns
  traits_layer <- fish_traits[species_in_layer, , drop = FALSE] %>%
    select(where(~ length(unique(.)) > 1))  
  
  # If no variable lines remain, we move on to the next layer
  if (ncol(traits_layer) == 0) next  
  
  # Correct formatting of biomass_layer with depth layer in rownames
  # Converts biomass data to matrix, which is required by FD::dbFD()
  biomass_layer <- matrix(as.numeric(depth_fish_biomass[layer, species_in_layer, drop = FALSE]), 
                          nrow = 1, dimnames = list(layer, species_in_layer))
  
  # Calculation of observed indices with FD::dbFD
  dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,
                          corr = corr_method)
  
  # Storage of observed values
  dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic
  dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis
  dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve
  dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv  
  
  # Simulations for each layer
  for (sim in 1:n_simulations) {
    randomized_traits <- randomize_traits(traits_layer)
    
    dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,
                                corr = corr_method)
    
    # Storage of simulated values
    dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic
    dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis
    dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve
    dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv 
  }  
}

# Calculation of the SES (Standardized Effect Size) ----
SES_results <- data.frame(depth_layer = rownames(dbFD_result_obs))

for (index in indices) {
  meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)
  sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)
  SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD
}

# Plot ----
results_df <- tidyr::pivot_longer(
  SES_results,
  cols = -depth_layer,
  names_to = "index",
  values_to = "SES_FD"
)

results_df$depth_layer <- factor(
  results_df$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
)

results_df$index <- factor(
  results_df$index,
  levels = c("FRic",
             "FDis",
             "FDiv",
             "FEve"),
  labels = c(
    "Functional richness",
    "Functional dispersion",
    "Functional divergence",
    "Functional evenness"
  )
)

ggplot(results_df, aes(x = depth_layer, y = SES_FD, color = depth_layer)) +
  geom_point(size = 3) +
  facet_wrap(~index) +
  geom_hline(yintercept = c(1.96, -1.96), linetype = "dashed", color = "gray40", size=0.8) +
  scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  labs(x = "", y = "Standard Effect Size (SES)") +
  theme_light() +
  theme(axis.text.x = element_blank(),
        strip.text.x = element_text(size = 14, color = "black"),
        strip.background = element_rect(fill = "white"),
        axis.title = element_text(size = 13),
        axis.text = element_text(size = 13)) +
  guides(col="none", fill="none")

ggsave("SES_dbFD.png", path = "figures", dpi = 700, height = 5, width = 7)

3.2. Normality and symmetry tests

• (Botta-Dukát, 2018)

• On the diversity index values obtained at random

# normality and symetry tests ----
# test_norm  <- SES_values_reps %>%
#   tibble::remove_rownames() %>%
#   filter(index %in% c("fdis", "fdiv", "feve", "fric"))
# 
# # Compute skewness and Lilliefors tests to check normaltiy and symetry 
# stat_indices <- test_norm %>% 
#   group_by(index, depth_layer) %>% 
#   summarise(
#     skewness = round(e1071::skewness(SES),2),
#     Lilliefors_pvalue = round(nortest::lillie.test(SES)$p.value,2)
#   )
# 
# htmltools::tagList(DT::datatable(stat_indices))

3.4. Two levels (stations and depth layer)

# # Model parameters ----
# n_simulations <- 9
# 
# # Correction method for non-Euclidean distances
# corr_method <- "lingoes" 
# depth_layers <- rownames(station_sp)  # Depth layers
# 
# # Indices to calculate
# indices <- c("FRic", "FDis", "FDiv", "FEve")  
# 
# # Data definition: list with depth and station matrices ----
# analysis_levels <- list(
#   depth_layer = depth_fish_biomass,
#   station = station_sp
# )
# 
# # Store results ----
# all_results <- list()
# 
# for (level in names(analysis_levels)) {
#   data_matrix <- analysis_levels[[level]]
#   layers <- rownames(data_matrix)
#   
#   # Matrices for observed and simulated results ----
#   dbFD_result_obs <- matrix(NA, nrow = length(layers), ncol = length(indices), 
#                             dimnames = list(layers, indices))
#   
#   dbFD_result_sim <- array(NA, dim = c(length(layers), n_simulations, length(indices)), 
#                            dimnames = list(layers, paste0("Sim.", 1:n_simulations), indices))
#   
#   # Loop through each layer ----
#   for (layer in layers) {
#     species_in_layer <- colnames(data_matrix)[data_matrix[layer, ] > 0]
#     traits_layer <- fish_traits[species_in_layer, , drop = FALSE] %>%
#       select(where(~ length(unique(.)) > 1)) %>%
#       droplevels()
#     
#     if (ncol(traits_layer) == 0) next
#     
#     biomass_layer <- matrix(as.numeric(data_matrix[layer, species_in_layer, drop = FALSE]), 
#                             nrow = 1, dimnames = list(layer, species_in_layer))
#     
#     # Calculate FD indices for the observed data ----
#     dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer, 
#                             messages = FALSE, corr = corr_method)
#     
#     # Store observed results
#     dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic
#     dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis
#     dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve
#     dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv  
#     
#     # Perform random simulations ----
#     for (sim in 1:n_simulations) {
#       randomized_traits <- randomize_traits(traits_layer)
#       
#       dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer, 
#                                   messages = FALSE, corr = corr_method)
#       
#       # Store simulated results
#       dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic
#       dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis
#       dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve
#       dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv
#     }
#   }
#   
#   # Calculate SES (Standardized Effect Size) ----
#   SES_results <- data.frame(level = level, name = rownames(dbFD_result_obs))
#   for (index in indices) {
#     meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)
#     sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)
#     SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD
#   }
#   
#   all_results[[level]] <- SES_results
# }
# 
# # Final result storage ----
# SES_results_depth_layer <- all_results$depth_layer %>% 
#  tidyr::pivot_longer(cols = -c(name, level), names_to = "index", values_to = "SES_FD") %>% 
#   mutate(depth_layer=name)
# 
# SES_results_station <- all_results$station %>%
#   rename(station=name) %>% 
#   tidyr::pivot_longer(cols = -c(station, level), names_to = "index", values_to = "SES_FD") %>% 
#   inner_join(metadata %>% select(station, depth), by = "station") %>%
#   # add column with depth layer
#   mutate(
#     depth_layer = case_when(
#       between(depth, 0, 174) ~ "Epipelagic",
#       between(depth, 175, 699) ~ "Upper mesopelagic",
#       between(depth, 700, 999) ~ "Lower mesopelagic",
#       between(depth, 1000, 2000) ~ "Bathypelagic"
#     )
#   ) 
# 
# # Combine datasets ----
# SES_results_combined <- bind_rows(SES_results_station, SES_results_depth_layer)
# 
# SES_results_combined$depth_layer <- factor(
#   SES_results_combined$depth_layer,
#   levels = c(
#     "Epipelagic",
#     "Upper mesopelagic",
#     "Lower mesopelagic",
#     "Bathypelagic"
#   )
# )
# 
# SES_results_combined$index <- factor(
#   SES_results_combined$index,
#   levels = c("FRic", "FDis", "FDiv", "FEve"),
#   labels = c(
#     "Functional richness",
#     "Functional dispersion",
#     "Functional divergence",
#     "Functional evenness"
#   )
# )
# 
# # Plot ----
# ggplot(SES_results_combined, aes(x = depth_layer, y = SES_FD)) +
#   geom_boxplot(data = SES_results_combined %>% filter(level == "station"),
#                aes(color = depth_layer, fill = depth_layer),
#                alpha = 0.08, width = 0.5, outlier.shape = NA) +
#   geom_point(data = SES_results_combined %>% filter(level == "depth_layer"),
#              aes(color = depth_layer),
#              size = 5, shape = 18) +
#   geom_jitter(data = SES_results_combined %>% filter(level == "station"),
#               aes(color = depth_layer),
#               size = 1, width = 0.2, alpha = 0.35) +
#   geom_hline(yintercept = c(1.96, -1.96), linetype = "dashed", color = "gray40", size = 0.8) +
#   scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
#   scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
#   facet_wrap(~ index, labeller = labeller(index = c(
#     FRic = "Functional richness", 
#     FDis = "Functional dispersion", 
#     FDiv = "Functional divergence", 
#     FEve = "Functional evenness"
#   ))) + 
#   labs(x = NULL, y = "Standard Effect Size (SES)") + 
#   theme_light() +
#   theme(
#     axis.text.x = element_blank(),  
#     strip.text.x = element_text(size = 14, color = "black"), 
#     strip.background = element_rect(fill = "white"),  
#     axis.title = element_text(size = 13),
#     axis.text = element_text(size = 13)
#   ) +
#   guides(col = "none", fill = "none")
# #ggsave("SES_dbFD_two_levels.png", path = "figures", dpi = 700, height = 7, width = 9)

• the red dashed line = normal distribution

# test_norm$depth_layer <- factor(
#   test_norm$depth_layer,
#   levels = c(
#     "Epipelagic",
#     "Upper mesopelagic",
#     "Lower mesopelagic",
#     "Bathypelagic"
#   )
# )
# 
# # Density plot with normal density overlay
# ggpubr::ggdensity(test_norm, x = "SES", fill = "gray", col = "white") +
#   facet_grid(depth_layer~index, scales = "free") +
#   ggpubr::stat_overlay_normal_density(color = "red", linetype = "dashed", size = 1) +  
#   theme_light() +
#   labs(title = "Distribution of values of the random communites in compariaosn to the normal distribution",
#        x = "Values", 
#        y = "Density")

4. CWM

4.1. CWM booststrap by depth layer

• each trawl is a replica

• non parametric Bootstrap

# Trait 
trait_boot <- morpho_data%>% 
  inner_join(metadata) %>% 
  select(-c(individual_code, years, longitude_start,
            latitude_start, longitude_end, longitude_end,
            volume_filtered, distance_trawled_m)) %>% 
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    depth,
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume,
    ventral_photophores, 
    gland_head,
    chin_barbel, 
    small_teeth, 
    large_teeth, 
    fang_teeth, 
    retractable_teeth, 
    internal_teeth
  ) %>%
  mutate_at(vars(ventral_photophores, 
                 gland_head,
                 chin_barbel, 
                 small_teeth, 
                 large_teeth, 
                 fang_teeth, 
                 retractable_teeth, 
                 internal_teeth), 
            funs(ifelse(. == "P", 1, ifelse(. == "A", 0, .)))) %>% 
  mutate(across(all_of(c("ventral_photophores", 
                         "gland_head",
                         "chin_barbel", 
                         "small_teeth", 
                         "large_teeth", 
                         "fang_teeth", 
                         "retractable_teeth", 
                         "internal_teeth")), as.numeric)) %>% 
  tidyr::pivot_longer(!c(species,depth), names_to = "trait", values_to = "values")%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

# Community 
community <-  rbind(data_biomass_2002_2019, data_biomass_2021_2022)%>%
  as.data.frame()%>%
  left_join(metadata) %>%
  select(species, biomass_sp, depth, volume_filtered)%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_sp))%>%
  select(-c(biomass_sp))%>%
  distinct()%>%
  select(-c(volume_filtered)) %>% 
  filter(biomass>0)

trait_filling <- traitstrap::trait_fill(
  # input data (mandatory)
  comm = community,
  traits = trait_boot,
  
  # specifies columns in your data (mandatory)
  abundance_col = "biomass",
  taxon_col = "species",
  trait_col = "trait",
  value_col = "values",
  
  # specifies sampling hierarchy
  scale_hierarchy = c("depth_layer", "depth"),
  
  # min number of samples
  min_n_in_sample = 5
)

# run nonparametric bootstrapping
np_bootstrapped_moments <- traitstrap::trait_np_bootstrap(
  trait_filling, 
  nrep = 100
)

sum_boot_moment <- trait_summarise_boot_moments(
  np_bootstrapped_moments
) %>% 
  mutate(trait= gsub("_"," ", trait)) %>% 
   mutate(trait=stringr::str_to_sentence(trait)) 

sum_boot_moment$depth_layer <- factor(
  sum_boot_moment$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
) 

# order traits 
sum_boot_moment$trait <- factor(
  sum_boot_moment$trait,
  levels = c(
    "Caudal throttle width",
    "Oral gape surface",
    "Gill outflow",
    "Large teeth",
    "Eye size",
    "Orbital length",
    "Small teeth",
    "Transversal shape",
    "Body depth",
    "Dorsal fin insertion",
    "Eye position",
    "Oral gape shape",
    "Oral gape position",
    "Internal teeth",
    "Lower jaw length",
    "Pectoral fin position",
    "Ventral photophores",
    "Operculum volume",
    "Pectoral fin insertion",
    "Head length",
    "Chin barbel",
    "Fang teeth",
    "Gland head",
    "Retractable teeth"
  )
) 

ggplot(sum_boot_moment, aes(x = depth_layer, y = mean)) + 
  geom_point(alpha = 0.5, size = 1, position = position_jitter(width = 0.2), aes(col=depth_layer)) +
  geom_boxplot(aes(group=depth_layer, col=depth_layer, fill=depth_layer), alpha=0.1)+
  scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  facet_wrap(~trait, scales = "free", ncol = 4) +
  theme_light() +
  labs(y="Community biomass-weighted mean (CWM) ")+
  guides(col="none", fill="none")+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.grid.major = element_blank(),
        axis.title.y.left = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        axis.text.y = element_text(size = 15),
        strip.background = element_rect(fill = "white"),
        strip.text.x = element_text(size = 16, color = "black"))

ggsave("CWM_boot_sum.png", path = "figures", dpi = 800, height = 13, width = 12)

4.2. PCA CWM

sum_boot_moment_pca <- trait_summarise_boot_moments(
  np_bootstrapped_moments
) %>% 
  ungroup() %>% 
  mutate(trait= gsub("_"," ", trait)) %>% 
  mutate(trait=stringr::str_to_sentence(trait)) %>% 
  select(c(trait, depth_layer, mean)) %>% 
  group_by(
    depth_layer, trait
  ) %>% 
  summarise(median_value= median(mean)) %>% 
  distinct() %>% 
  ungroup() %>% 
  tidyr::pivot_wider(id_cols=depth_layer, values_from = "median_value", names_from = "trait") %>% 
  tibble::column_to_rownames(var = "depth_layer")


res.pca <- FactoMineR::PCA(sum_boot_moment_pca, graph = FALSE)

pca_var <- factoextra::fviz_pca_var(
  res.pca,
  repel = TRUE,
  pointsize = 2,
  arrowsize = 0.1,
  label = "var",
  title = "",
  col.circle = "gray80",
  col.var = "black",
  labelsize=5,
  alpha.arrow=0.5
) +
  labs(x = "", y = "") +
  theme_light() +
  theme(
    aspect.ratio = 1,
    panel.grid.minor =  element_blank(),
    panel.grid.major = element_blank())

pca_ind <- factoextra::fviz_pca_ind(res.pca, label="none", title="") +
  theme_light() +
  xlim(c(-7,7))+
  labs(x="", y="")+
  theme(aspect.ratio = 1)

combined_plot <- gridExtra::grid.arrange(pca_var, pca_ind, ncol = 2)

ggsave("PCA_CWM.png", plot = combined_plot, path = "figures", dpi = 700, height = 6, width = 11)

5. Functional rarity

• Uniqueness - geographical restrectiveness

which traits are mainly driving uniqueness ?

df <- fish_traits %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "species") %>%
  left_join(sp_ui) %>%
  tibble::column_to_rownames(var = "species")

# Identify numerical and categorical traits
numeric_traits <- colnames(df)[sapply(df, is.numeric) &
                                 colnames(df) != "Ui"]

categorical_traits <- colnames(df)[!sapply(df, is.numeric) &
                                     colnames(df) != "Ui"]

# Initialize a data frame to store all results
combined_results_df <- data.frame(
  Trait = character(0),
  Type = character(0),
  Eta_R_squared = numeric(0),
  P_value = numeric(0),
  stringsAsFactors = FALSE
)

# Step 1: Perform Kruskal-Wallis test for categorical traits
for (categorical_trait in categorical_traits) {
  kruskal_result <- kruskal.test(df$Ui ~ df[[categorical_trait]])
  
  # Calculate eta-squared for Kruskal-Wallis
  n_groups <- length(unique(df[[categorical_trait]]))
  n_total <- length(df$Ui)
  h_value <- kruskal_result$statistic
  eta_squared <- (h_value - (n_groups - 1)) / (n_total - n_groups)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = categorical_trait,
      Type = "Categorical",
      Eta_R_squared = as.numeric(format(eta_squared, scientific = FALSE)),
      P_value = as.numeric(format(kruskal_result$p.value, scientific = FALSE)),
      stringsAsFactors = FALSE
    )
  )
}

# Step 2: Fit linear models for numerical traits
for (numeric_trait in numeric_traits) {
  lm_result <- lm(Ui ~ df[[numeric_trait]], data = df)
  summary_stats <- summary(lm_result)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = numeric_trait,
      Type = "Numeric",
      Eta_R_squared = as.numeric(format(summary_stats$r.squared, scientific = FALSE)),
      P_value = as.numeric(format(
        summary_stats$coefficients[2, 4], scientific = FALSE
      )),
      stringsAsFactors = FALSE
    )
  )
}

# Round the numeric columns to three decimal places
combined_results_df[, c("Eta_R_squared", "P_value")] <-
  round(combined_results_df[, c("Eta_R_squared", "P_value")], 3)

htmltools::tagList(DT::datatable(combined_results_df))

6. Appendices

# Convert sp_faxes_coord_fish to a data frame
sp_coord_community_2 <- as.data.frame(sp_faxes_coord_fish[, c("PC3", "PC4")]) %>%
  tibble::rownames_to_column(var = "species")

# Create a list to store data frames for each depth layer
all_data_layers_2 <- lapply(unique(sp_all_layers$depth_layer), function(layer) {
  sp_layer <- sp_all_layers %>%
    filter(depth_layer == layer)
  
  #presence or absence
  sp_coord_layer <- sp_coord_community_2 %>%
    mutate(layer_presence = ifelse(species %in% sp_layer$species, "yes", "no")) %>%
    mutate(depth_layer = layer)
  
  return(sp_coord_layer)
})

# Combine all layers into one data frame
sp_coord_community_all_layers <- bind_rows(all_data_layers_2)

# Function to calculate the convex hull for a given data frame
calculate_hull <- function(data) {
  if (nrow(data) < 3) return(data)  
  data %>%
    slice(chull(PC3, PC4))
}

# Initialize lists to store hulls
hull_all_combined <- NULL
hull_layer_combined <- NULL

for(layer in unique(sp_coord_community_all_layers$depth_layer)) {
  sp_layer <- sp_coord_community_all_layers %>%
    filter(depth_layer == layer)
  
  hull_all <- calculate_hull(sp_layer) %>%
    mutate(depth_layer = layer, type = "all")
  
  hull_layer <- calculate_hull(sp_layer %>% filter(layer_presence == "yes")) %>%
    mutate(depth_layer = layer, type = "present")
  
  hull_all_combined <- bind_rows(hull_all_combined, hull_all)
  hull_layer_combined <- bind_rows(hull_layer_combined, hull_layer)
}

# Define depth layer levels
depth_levels <- c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic")

sp_coord_community_all_layers$depth_layer <- factor(sp_coord_community_all_layers$depth_layer, levels = depth_levels)
hull_all_combined$depth_layer <- factor(hull_all_combined$depth_layer, levels = depth_levels)
hull_layer_combined$depth_layer <- factor(hull_layer_combined$depth_layer, levels = depth_levels)

# Define colors for each depth layer
depth_colors <- c("Epipelagic" = "#FEA520", 
                  "Upper mesopelagic" = "#D62246", 
                  "Lower mesopelagic" = "#6255B4", 
                  "Bathypelagic" = "#3C685A")

# Plot
plot_depth <- ggplot() +
  # Plot present species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "yes"), 
             aes(PC3, PC4, col = depth_layer), size = 2, shape=19) +
  # Plot absent species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "no"), 
             aes(PC3, PC4), col = "grey", shape = 8, size = 2) +
  # Plot hulls
  geom_polygon(data = hull_all_combined, aes(x = PC3, y = PC4), 
               fill = "gray70", alpha = 0.2) +
  geom_polygon(data = hull_layer_combined, aes(x = PC3, y = PC4, group = depth_layer, fill = depth_layer), 
               alpha = 0.4) +
  theme_light() +
  scale_color_manual(values = depth_colors) +
  scale_fill_manual(values = depth_colors) +
  guides(col = "none", shape = "none", fill = "none") +
  facet_wrap(~ depth_layer, nrow=1 ) +  
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )


beta_fd_indices_fish <- mFD::beta.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC3", "PC4", "PC3", "PC4")],
  asb_sp_occ       = asb_sp_fish_occ,
  check_input      = TRUE,
  beta_family      = c("Jaccard"),
  details_returned = T)

Serial computing of convex hulls shaping assemblages with conv1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of convex hulls shaping assemblages with conv2

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_geom_coord

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_rcdd_coord & qhull.opt1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_rcdd_coord & qhull.opt2

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%

vertices_community <-beta_fd_indices_fish$"details"$"asb_vertices"$Bathypelagic 

sp_coord_community_2 <- as.data.frame(sp_faxes_coord_fish[, 3:4]) %>% 
  tibble::rownames_to_column(var = "species") %>% 
  mutate(vertices=case_when(species%in%vertices_community~"vertice",
                            !species%in%vertices_community~"not_vertice")) 

hull <- sp_coord_community_2%>% 
  slice(chull(PC3, PC4))

plot_com <- ggplot(data = sp_coord_community_2,aes(PC3, PC4)) +
  scale_color_manual(values = c("grey", "black"))+
  geom_point(size = 2.5,aes( shape=vertices, col= vertices)) +
  geom_polygon(data=hull, aes(x = PC3, y = PC4), 
               fill = "gray70", alpha = 0.2) +
  theme_light() +
  geom_hline(yintercept = 0, linetype = 2, col = "gray45") +
  geom_vline(xintercept = 0, linetype = 2, col = "gray45") +
  guides(col = "none", shape = "none", fill = "none") +
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )

ggpubr::ggarrange(
  labels = c("A","B"),
  plot_com,   
  plot_depth, 
  nrow = 2
) 

ggsave("functional_space_PC3_4.png", path = "figures", dpi = 800, height = 8, width = 10)